NASA Contract Report 
NAS 1 -00086-A003 



Computation of Sound Propagation by Boundary Element 
Method 


Yueping Guo 

The Boeing Company 
Mail Code H013-B308 
5301 Bolsa Avenue 
Huntington Beach, CA 92647 

vueping.guo@boeing.com 


September 2005 


National Aeronautics and Space Administration 
Langley Research Center 
Hampton, VA 23681 



Table of Contents 


1. Introduction 5 

2. Governing Equations 7 

3. Numerical Scheme 11 

4. Sound Scattering by a Sphere 14 

5. Sound Reflection by a Plate with Uniform Flow 17 

6. Sound Propagation over a Hump 20 

7. Far Field Formulation 22 

8. Calculation of Derivatives 23 

9. Conclusions 25 

Nomenclature 26 

Fist of Figures 27 

References 28 


2 




Acknowledgments 


The work reported here was conducted under NASA Contract NAS 1-00086, under the Quiet 
Aircraft Technology (QAT) program. The author would like to thank the task monitor, Dr. 
Mehdi Khorrami of NASA LaRC, for his support and encouragement. 


3 




Abstract 


This report documents the development of a Boundary Element Method (BEM) code for 
the computation of sound propagation in uniform mean flows. The basic formulation and 
implementation follow the standard BEM methodology; the convective wave equation and the 
boundary conditions on the surfaces of the bodies in the flow are formulated into an integral 
equation and the method of collocation is used to discretize this equation into a matrix equation 
to be solved numerically. New features discussed here include the formulation of the additional 
terms due to the effects of the mean flow and the treatment of the numerical singularities in the 
implementation by the method of collocation. The effects of mean flows introduce terms in the 
integral equation that contain the gradients of the unknown, which is undesirable if the gradients 
are treated as additional unknowns, greatly increasing the sizes of the matrix equation, or if 
numerical differentiation is used to approximate the gradients, introducing numerical error in the 
computation. It is shown that these terms can be reformulated in terms of the unknown itself, 
making the integral equation very similar to the case without mean flows and simple for 
numerical implementation. To avoid asymptotic analysis in the treatment of numerical 
singularities in the method of collocation, as is conventionally done, we perform the surface 
integrations in the integral equation by using sub-triangles so that the field point never coincide 
with the evaluation points on the surfaces. This simplifies the formulation and greatly facilitates 
the implementation. To validate the method and the code, three canonic problems are studied. 
They are respectively the sound scattering by a sphere, the sound reflection by a plate in uniform 
mean flows and the sound propagation over a hump of irregular shape in uniform flows. The first 
two have analytical solutions and the third is solved by the method of Computational 
Aeroacoustics (CAA), all of which are used to compare the BEM solutions. The comparisons 
show very good agreements and validate the accuracy of the BEM approach implemented here. 
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1. Introduction 

In practical applications such as airframe noise prediction, noise propagation very often 
involves complex geometry, the effects of which can only be dealt with by numerical methods. 
There are various numerical methods studied in the past, including the Boundary Element 
Method (BEM). Though the basic theories of these methods are well established, their 
applications to practical problems are limited by their heavy requirements on computation 
resources. Compared with other numerical methods, BEM has the advantageous feature of only 
involving surface discretization, making the numerical computation less demanding, both in 
computation time and in data storage. Because of this, BEM plays an important and useful role in 
practical applications of noise prediction. In this report, we document the formulation, 
implementation and validation of a BEM code. 

The Boundary Element Method has been an active research topic for many years and the 
fundamental theory and formulation have been well established and documented in the open 
literature (e.g. Ref 1 to Ref 7). Thus, it is appropriate to point out that the work reported here is 
not intended to provide any new break through in this topic. The objective of this report is 
twofold. The first is to provide a document for the implementation and validation of this method 
with all the relevant technical details so that it serves as a user manual. Our second objective is to 
discuss a few features related to the implementation of BEM that do not seem to have reported in 
the open literature. These include a treatment of the mean flow effects to eliminate the gradients 
of the unknown in the integral equation, a discretization scheme that is free of numerical 
singularity and a formulation that computes the high order spatial derivatives of the solution. 
Again, these features are not fundamental break through of BEM, but they facilitate its 
implementation and provide incremental improvement to the method. 

The formulation of BEM is well established. It involves converting the governing 
equation, the convective wave equation for applications with uniform mean flow, and the 
boundary condition, the vanishing of the normal derivative of the unknown on the surfaces of 
bodies embedded in the flow, into an integral equation on the body surfaces through the 
application of the Greens theorem. The integral equation is solved numerically by discretizing it 
into a matrix equation, which leads to solutions of the unknown on the surfaces. Solutions in the 
field can then be found from these surface solutions by integration. In the integral equation 
formulation, both the unknown and its gradients appear in the integrand in cases of non-zero 
uniform mean flow. Traditionally, this is dealt with by one of the two approaches. One is to 
differentiate the integral equation to derive an additional set of equations so that both the 
unknown and its gradients are solved simultaneously by two sets of coupled matrix equations. 
Clearly, compared with the cases without mean flow where only the unknown itself is solved, 
this approach significantly increases the computation time and storage requirement, by doubling 
the dimension, and hence quadrupling the size, of the matrix equation to be solved. The second 
approach to deal with the gradients of the unknown in the integral equation is to approximate the 
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gradients numerically, by finite differencing, for example, so that the integral equation reduces to 
one that only contains the unknown itself. While this approach does not increase the dimensions 
of the matrix equation to be solved, it relates the unknown at a surface point to those at its 
neighboring points in the integral equation, making the discretization scheme and coding more 
complex. It also introduces additional numerical error in approximating the derivatives by 
numerical differencing. 

In this report, we will present a different approach that avoids the disadvantages of the 
two previous approaches. We will show that the terms involving the gradients of the unknown in 
the integral equation can be reformulated in terms of the unknown itself. By this reformulation, 
the integral equation to be solved only involves the unknown itself and the structure of the 
integral equation is very similar to that for cases without mean flow. The reformulation is exact 
and does not involve any numerical differentiation. This is clearly advantageous and desirable, 
both ensuring the accuracy of the approach and minimizing the requirement of computation 
resources. 

There are many ways to discretize the integral equation in BEM for numerical solution, 
as extensively discussed and review in the literature. The methods range from the collocation 
method to the spectral method, the former corresponding mathematically to the use of a local 
weight in the kernel of the integral equation and the latter to the use of a global weight. Each of 
these methods has its advantages and disadvantages. We choose to work with the collocation 
method because of its ease of implementation. To avoid its disadvantage of having to deal with 
numerical singularities when the solution point coincides with one of the surface integration 
points, we carry out the surface integration by further dividing each of the discretized individual 
surface elements into sub-triangles. When numerically calculating the surface integration on an 
element, by summing the products of the integrand values and the sub-triangle surface areas, the 
solution point and the surface integration points never coincide because the former is one of the 
vertexes of the triangles and the latter is their geometric centers. This eliminates the need for 
asymptotic analysis of the integral equation, as necessary in typical implementation of BEM by 
collocation method. 

One application of BEM is as a Greens function solver to compute the propagation 
effects of sound waves in the prediction of noise from unsteady flows. In such applications, it is 
often required to compute the derivatives of the Greens function (Ref 8 to Ref 10). It can in 
theory be done by numerical differentiation once the Greens function itself is found in all spatial 
domains. This, however, is not the preferred approach, both because of the extra computation 
time required to carry out the numerical differentiation and because of the accumulated 
numerical errors in the final results. The latter is probably the more serve limitation of the two, 
especially with non-uniform spatial grid and low order schemes of numerical differentiation. 
These limitations can be overcome by computing the high order derivatives of the Greens 
function directly from the solutions of the Greens function on the body surfaces, namely, the 
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solutions of the BEM matrix equation. In this report, we derive the formulas for computing the 
high order derivatives, which, though not difficult, is quite involved and tedious. 

To validate the BEM approach and the implemented code, we choose to study three 
canonical problems. The first is the sound scattering from a rigid sphere in a static medium. 
Clearly, this problem is considered here because of the availability of analytical solutions, as 
studied and documented extensively in the literature. The analytical solutions provide a good 
validation with minimum uncertainty. The problem, however, is only for static medium without 
mean flows. For problems with both mean flows and boundary surfaces present, exact analytical 
solutions basically do not exist. Our validation will then rely on approximate solutions and on 
comparisons with other numerical methods. For this purpose, we also use a numerical code, 
based on the methods of Computational Aeroacoustics (CAA), to study the second and third 
canonic problem, which are respectively the sound reflection from a large flat plate, placed in 
parallel with the direction of a uniform mean flow, and the sound propagation over a hump in a 
uniform mean flow. For a plate whose dimensions are large compared with the sound 
wavelength, the solutions away from the plate edges are not significantly affected by the edges 
so that the analytical solution for an infinity plate serves as a good approximation to the finite 
plate solution. By studying this problem by both CAA and BEM, and comparing the results with 
the approximate analytical solutions, not only the BEM code gets validated for a case with mean 
flow, but also the relevance of the CAA code is established so that it can provide comparisons in 
cases where even approximate solutions are not available. Our third validation is such a case, 
where sound waves propagate over an irregularly shaped hump in a mean flow. In this case, the 
validation of BEM solutions entirely replies on the comparisons between BEM and CAA. As 
will be seen in subsequent sections, all three validation cases show good accuracy of the BEM 
solutions. 

2. Governing Equations 

When applying BEM to acoustic problems, we seek to solve the convective wave 
equation in frequency domain in the form of, 


where g is the unknown to be solved, which can be regarded as the velocity potential, and q can 
be any source distribution. We have introduced k (} and M to respectively denote the acoustic 
wavenumber and the vector mean flow Mach number. The acoustic wavenumber k o is related to 
the angular frequency co and sound speed c 0 by 


(-ik 0 +M-V) 2 g-V 2 g = 4(x), 


( 2 . 1 ) 


k 0 =co/c Q , 


(2.2) 


and the vector Mach number M is defined by the flow velocity U by 
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(2.3) 


M = U/c 0 . 

The governing equation (2.1) is defined by the Cartesian coordinate system 

x = {x, , x 2 , x 3 }, (2.4) 

which is fixed on the rigid bodies embedded in the flow. The geometry is illustrated in Figure 1. 
On the boundary surfaces of the bodies, the boundary condition is the vanishing of the normal 
derivative of g, namely, 


dn 


= n ■ Vg = 0, 


(2.5) 


where n denotes the unit normal of the surfaces, 
pointing into the flow, as illustrated in Figure 1. 

From this boundary condition, it is also clear that 
the unknown g is the velocity potential and its 
gradient is the velocity vector. This unknown 
quantity is denoted here by g because the 
solutions of the governing equation (2.1), with 
the boundary condition (2.5), become the Greens 
function if the source distribution is a 
concentrated point source, which is one main 
objective in implementing the BEM here. 

The standard procedure of using BEM to 
solve the governing equation (2.1) is to apply the 
Greens theorem to derive an integral representation of the solution, which starts with defining an 
adjoint function go by 

(- \k 0 - M • V) 2 g 0 - V 2 g 0 = S(x - z). (2.6) 

Apparently, this equation differs from the governing equation (2.1) with a reversed flow Mach 
number and with the replacement of the source term by the delta function on the right hand side, 
which defines a point source located at z. The next step is to multiply (2.1) by go and (2.6) by g, 
and to take the difference of these two results, which leads to 

V ■ {- 2ik 0 Mgg 0 + M(g 0 Vg - gVg 0 ) ■ M - g 0 Vg + gVg 0 } 

= g 0 (x-z)g(x)-g(x)£(x-z). 

By integrating this over the flow domain outside the embedded bodies, it is straightforward to 
derive 



Figure 1 Illustration of the geometry and 
coordinate system. 


J[2i k 0 M n g 0 g + (n -M n M) ■ (g 0 Vg - gVg 0 )]ds = Jg 0 (z - x)q(z)dv - Cg(x). (2.8) 

S(z) V(z) 
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Here the left hand side follows from the use of the divergence theorem, which reduces the 
volume integration over the left hand side of (2.7) to a surface integration. The surfaces of the all 
the rigid bodies in the flow are collectively denoted by S(z), with differential surface element ds, 
and the component of the mean flow Mach number in the surface normal direction is denoted by 
M„, namely, 

M n = n M. (2.9) 

The right hand side of the result (2.8) contains a volume-integration over the source distribution 
q, with the differential volume element denoted by dv. The last term in the right hand side results 
from the trivial application of the delta functions to carry out the volume-integration. In carrying 
out this integration, the coefficient C is introduced to account for the three cases where the field 
point x is respectively inside, on the boundary surface and outside the flow. Its values are given 
by 


fl 


C = 


0.5 

0 


for x in flow 

for x on boundary surface 
for x outside flow 


(2.10) 


which is well-known in BEM formulations and accounts for the singular nature of the surface 
integral when the field points are on the surface. 

Now the boundary condition (2.5) can be applied to simplify the result (2.8) by 
eliminating terms involving the normal derivatives of g, which leads to 

Cg(x)= jg 0 (z-x)q(z)dv 

V(z) (2.11) 

- j [g(2ik 0 M„g 0 + M „ M ■ Vg 0 - n -Vg 0 )~ g 0 M n M ■ Vg]ds. 

S(z) 

It should be noted that the gradients inside the surface integration are with respect to the 
coordinate z. This is the integral equation to be solved for the function g. It can be seen that the 
integral on the right hand side not only contains the unknown function g, but also its gradients, 
respectively in the two terms in the integrand. Apparently, the gradient term is due to the 
presence of the uniform mean flow since the term involving the gradient of the unknown is 
proportional to the mean flow Mach number. In a static acoustic medium, this term vanishes and 
the integral is entirely given by the function g itself. 

To solve the equation (2.11) numerically, it is desirable to have only g as the unknown 
without the term involving its derivatives. This term can be further manipulated by writing the 
vector Mach number M and the gradient of g in a body-fitted coordinate system that has two 
components on the body surface and the surface normal as its third component. In this coordinate 
system, we have 
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( 2 . 12 ) 


M-Vg=M n ^ + M s -V s g 
an 


M s -V s g, 


where the last step follows from the use of the boundary condition (2.5) and the subscript s 
denotes the two coordinates in the two tangential directions on the surface ( s = 1, 2). These 
surface coordinates can be defined as the arc lengths on the surface, but in vector form, a surface 
quantity can be simply written as the difference between the full vector and its normal 
component. Thus, we have the definitions 


M ( = M - A/ n n and 



From (2.12), the last term in (2.11) becomes 


(2.13) 


J So M „ M • Vgds = J g 0 M„M s -V S gds 

S(z) S(z) 

= f[ v s •(^o M „M s g)-gV s -( go M n M s )]ds. 

S( z) 


(2.14) 


The first term in the above result is in a divergence form, in terms of the surface coordinates, 
which integrates to zero for finite bodies with closed surfaces. The last term only involves g as 
the unknown. Thus, the integral equation (2.11) reduces to the compact form 


Cg(x) + f A(z,x)g(z)ds = <2(x), 

S( z) 


(2.15) 


where Q is written for 

Q= \g 0 (z-x)q(z)dv, (2.16) 

V (z) 

and the quantity A denotes 

A = 2« 0 M„ S o+(M„M-n).V So +V,.( So M„M,). (2.17) 

Apparently, the result (2.15) is an integral equation with g as the only unknown. 

To proceed from here, the adjoint function go must be found, as the solution to (2.6) in an 
infinity space. The solution to this equation can be found easily by standard procedures ( e.g . Ref 
1 1 and 12) in the form of 


where R denotes 


g 0 (z-x) = 


4/rR 4 


ii 0 (fi‘+M.(z-x))//? 2 


R = +((z-x)-M) 2 , 


(2.18) 


(2.19) 
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with P defined by 


P 2 =\-M 2 . (2.20) 

Here M is the mean flow Mach number. 

From the solution (2.18), the gradients of gO with respect to z can be found easily as 


VSc 


J i* 0 (R*+M (z-x ))/^ 2 I ■ K 


VR* 


AnR* 


f \{vR + m) 

p 1 1 J R 


Clearly, this can be symbolically denoted in the compact form 


(2.21) 


with the vector F denoting 


Vg 0 = So F > 


(2.22) 


F = i 4r{ VR * +M F^ _ - C- 23 ) 

P K 

The gradients of go with respect to z are needed in calculating the surface integral in the quantity 
defined by (2.17). 

3. Numerical Scheme 

To solve the integral equation (2.15), we divide the surface S into N surface elements, 
each having Kj sides (j = 1, 2, 3, ... AO- In theory, the elements can have any number of sides, but 
in practical applications, they are most likely to be triangles and quadrilaterals (Kj = 3 or 4). By 
dividing the body surfaces into elements, the surface integration in (2.15) is to be performed on 
each small element, on which, the unknown g( z) is assumed to be constant so that it can be 
moved outside the integration. This leads to 

Cg(x) + Y J g(^j)lM^^)ds = Q(x), (3.1) 

j= 1 Sj 

where z j is the coordinates of the geometric centers of the jth panel. The surface integral is now 
over the small elements and only involves the known quantity A, namely, 

J A(x,z)ds = j{2ik 0 M n g 0 +(M n M-n)-Vg 0 +V s -{g 0 M n M s )}ds. (3.2) 

Sj Sj 

The evaluation of this term is straight forward once the singularities in the free space Greens 
function go are dealt with, which is conventionally done by asymptotic analysis in the method of 
collocation. The asymptotic analysis can be avoided by suitably treating the integration on the 
elements, as described in the following. 
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To this end, we start with the first two terms in the integrand of (3.2) and divide the 
element Sj into Kj sub-triangles by connecting each vertex of the element with its geometric 
center zy, as illustrated in Figure 2. The integration over an individual sub-triangle can then be 
done by evaluating the integrand function at the geometric center of the triangle and multiplying 
the result by its area. The total integration over the jth element, namely, over Sj, is of course the 
summation over all the sub-triangles. For the first two terms in the integrand of (3.2), this leads 
to the result 


Y^{2ik 0 M n g 0 (x ~z k ) + (M n M -n) • Vg 0 (x - z k )}s t , (3.3) 

k= 1 

where the summation is over the number of sub- 
triangles in the element, which equals to the 
number of sides. For the kxh sub-triangle (k = 

1 ...Kj), its geometric center and surface area are 
respectively denoted by z k and s k . Both are 

illustrated in Figure 2. 

For the last term in the integrand of (3.2), 
the divergence can be used to reduce the surface 
integration to line integration, along the sides of 
the element. Each line integral can then be 
evaluated by the value of the integrand function 
at the center of the line segment multiplied by 
the length of the segment. This leads to 

K i 

Zso( x -z A .)M fI M 5 ■n k a k , (3.4) 

k = 1 

where z k is the center of the k\h side of the element, n /f is its normal, in the same plane as the 

sub-triangle and pointing outward away from the element, and a k is its length. Again, these 
definitions are illustrated in Figure 2. 

By collecting (3.3) and (3.4), the surface integration in (3.1) becomes a summation over 
the number of sides of the jth panel in the form of 

k j 

\Ads = ^[{2ifc 0 M fI g 0 (x-z it ) + (M fl M-n)-Vg 0 (x-z,)}^ + g 0 (x - z k )M n M s ■ n k a J (3.5) 

Sj k=l 

In this result, the quantities that depend on the individual sub-triangles are explicitly indicated by 
the subscript k. The quantity M n , defined by (2.9), is also evaluated for each sub-triangle because 
the surface normal can vary within each element. Thus, the surface normal n should be evaluated 
for each sub-triangle in the surface integration. 


jth Surface Element kth Sub-Triangle 



z k = Mi triangle center 
z k = kth side center 
n, = Mi side normal 

Figure 2 Geometry and definitions of a 
surface element and its sub-triangles 
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The result (3.5) gives the surface integration on the jth surface element, whose surface is 
denoted by Sj and on which the solution of g is assumed to be constant and represented by its 
value at the geometric center of the surface element zy. The result is valid for any point x so that, 
when substituted into (3.1), it can be used to compute field solutions as well as solutions on the 
surfaces. The matrix equation to solve for g can be constructed by setting (3.1) to the geometric 
centers of the surface elements. This leads to a matrix equation in the form of 

\gi + Yj a ijgj=Qi’ ( 3 - 6 ) 

1 j = i 

which unknowns of the equation are gj, defined by 

8i =£(X), (3.7) 

with Xj being the coordinates of the geometric centers of the z'th surface element. The source 
terms of (3.6) are similarly defined by 


Q i =Q(* i ), (3.8) 

and the coefficients of the matrix equation is given by 

a ij =jA(x ( .,z)ds, (3.9) 

s j 

namely, given by (3.5) with x replaced by x,-. It can be noted that when using (3.5) to calculate 
the coefficients of the matrix equation (3.6), the adjoint the free space Greens function has to be 
computed in the forms 


So( x ,~zJ and g 0 ( x i~ z k)- (3-10) 

Since the three position vectors in these respectively denote the geometric center of the z'th panel, 
the geometric center of the kth sub-triangle of the jth panel and the center of the kth side of the 
jth panel, the three points never coincide, even in the case of i = j, which is very desirable in 
numerical implementation because the computation of the matrix coefficients is always free of 
singularity. 

In comparison with various BEM formulations and implementations, the formulations 
given in this section and the previous section have two main attractive features. One is that the 
integral equation does not involve the gradients of the unknown, which limits the requirement on 
computation resources to the same as that for the cases without mean flow, both when solving 
the matrix equation for the solutions on the surfaces and when performing the surface integration 
to find solution in the field. The second attractive feature is that the coefficients of the matrix 
equation is always free of numerical singularities, even though the discretization of the integral 
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equation follows the approach of co-location, which usually requires asymptotic analysis to 
avoid singularities. These two features make the implementation very straightforward. 

4. Sound Scattering by a Sphere 

Sound scattering by a sphere is one of the 
very few model problems that have closed form 
solutions in a three-dimensional space with a finite 
size body imbedded in a static acoustic medium. 

Though the analytical solution does not include 
mean flow, it can serve as a good benchmark for 
numerical code validation. The solution is well 
documented in the literature and text books on 
acoustics (e.g. Ref 12). For a sphere of radius a 
with its center chosen to be the origin of a 
spherical coordinate system (r, 0, (f)), the sound 
scattering problem involving the sphere and a point 
source is axisymmetrical if the line between the 
sphere center and the point source is chosen to be 
the axis of 0= 0, as illustrated in Figure 3. Thus, the coordinate (j ) docs not appear in the solution 
and the location of the point source is specified by its radial distance r = r$. The governing 
equation to be solved is simply 



Figure 3 Geometry and coordinate system 
for sound scattering by a sphere. 


ko p + V 2 p = ~ — S(r - r o )S(0), 
r o 

subject to the boundary condition 


(4.1) 


dp 

dr 


= 0 


on 


r = a. 


(4.2) 


where p is the total acoustic pressure. The solution to this problem is well documented and can 
be written in the form 


P 


^ (2m + 1)-P (cos 6)\ <P - 
to ! 


j'n,( k 0 a ) 

h' m (k 0 a ) 


KAKr,)h(k,r). 


where <E> is used to denote 


\j m (k Q r 0 )h m (k Q r) for r > r 0 
\j m (k 0 r)h m (k 0 r 0 ) for r < r 0 ' 


(4.3) 


(4.4) 
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This term corresponds to the standard solution for 
a point source of angular frequency co, which 
determines the acoustic wavenumber by ko. The 
solution is in terms of the spherical coordinate 
system (r, 6, (/)). The expanded solution is 
expressed by the Lagrange polynomial P m of order 
m and the spherical Bessel and Hankel function of 
order m, respectively denoted by j m and h m . The 
prime over the spherical Bessel and Hankel 
functions indicates differentiation with respect to 
their argument. This solution can be easily 
computed with the summation truncated at some 
large number. The terms in the summation become 
progressively smaller so that the truncation only 
leads to a small error in the final results. 

The problem defined by (4.1) and (4.2) is 
solved by BEM to validate the BEM code (setting 
the mean flow to zero). Some comparisons are 
given in Figure 4 to Figure 7. In all the 
comparisons, the source location is set at ro = 4 a 
and the analytical solutions are computed from 
(4.3) and (4.4). Figure 4 shows comparisons at 
two field locations; the upper diagram is for a 
location in the insonified region on the same side 
of the sphere as the source at r = 5a and 8 = 0, and 
the lower diagram for a location in the shadow 
region at r = 5a and 6 = n. Both diagrams plot the 
real part of p (in red color), the imaginary part of p 
(in green) and its amplitude (in blue), as a function 
of the non-dimensional frequency kOa. The 
symbols are the results of the BEM computations 
and the curves are the analytical solutions. 
Apparently, the agreement between the two is very 
good, not only in the insonified region where the 
direct radiation from the source dominates, but 
also in the shadow region where the dominant 
noise is from the scattering of the sphere. This 
illustrates the accuracy of the BEM method 
implemented here. 





Figure 4 Comparisons between BEM 
(symbols) and analytical solutions (curves), 
with the upper and lower plot respectively 
for 6 = 0 and 8 = n, both for r = 5a. 
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Figure 5 Directivity comparisons for 
scattering by a sphere on the circle r = 1.5 a 
for k 0 a = 9.24. 
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Im(p): BEM 
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Re(p): BEM 
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Figure 7 Comparisons of the imaginary 
part of pressure between BEM and 
analytical solutions for k {) a = 9.24. 


Figure 6 Comparisons of the real part of 
pressure between BEM and analytical 
solutions for k 0 a = 9.24. 


To show global comparisons, Figure 5 and Figure 7 plot spatial distributions of the 
scattered pressure, both for a fixed k 0 a value of 9.24, which corresponds to a frequency of 500 
Hz for a sphere of radius of one meter. This relatively large value is chosen here to give 
sufficient variations on the sphere surface to test the capability of the BEM code in dealing with 
high frequency sources. Figure 5 plots the amplitude of p at field locations along a circle with r = 
1.5a, with the symbols for the results of BEM and the curve for the analytical solutions, which is 
commonly regarded as the directivity of the noise. Figure 6 and Figure 7 give a broad view 
comparison between the two solutions by plotting the real and the imaginary part of the pressure 
field, respectively on the surface of the sphere and the symmetry plane. Apparently, the two sets 
of solutions are almost identical. 


16 





5. Sound Reflection by a Plate with Uniform Flow 


When mean flows are present, closed 
form analytical solutions usually do not exist. 

Approximate analytical solutions may be found 
for simply geometries, one example being sound 
reflection by a large plate in a uniform flow. The 
solution is exact for an infinite plate, but for 
comparisons with numerical methods such as 
BEM and CAA which can only treat finite plates, 
the solution can be used in the region away from 
the edges of the plate. The sound field in this 
region is almost the same as that from an infinity 
plate and the analytical solution for an infinite 
plate provides a good approximation to the finite 
plate problem. 

The governing equation for the sound reflection problem is 

(-i k 0 +M-V) 2 p -V 2 p = S(x-x 0 )q(co), (5.1) 

where M is the mean flow Mach number and x 0 is the source location. Without loss of 
generality, we can choose the Cartesian coordinate system x such that two of the three axes are 
on the plate and the point source is on the axis perpendicular to the plate. The source location is 
then specified by the distance between the plate and the source, denoted by h, and the location 
vector becomes 

X 0 ={0,0,M. (5.2) 

Accordingly, the mean flow is in the positive x\ direction and the mean flow Mach number is 
given by 



Figure 8 Illustration of the geometry and 
coordinate system for the problem of sound 
reflection by a plate in a uniform flow. 


M = {M ,0, 0}. (5.3) 

In this coordinate system, the boundary condition is simply the vanishing of the normal 
derivative of p on the plate, expressed by, 

^- = 0 on * 3 = 0 . (5.4) 

OX 3 

The geometry and the coordinate system are illustrated in Figure 8 . 

In the governing equation (5.1), we have defined the source amplitude by a frequency 
dependent quantity, given by 
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(5.5) 


fl 


q(co) = 


cos- 

0 


n ( I k Q a I -1) 


I k Q a l< 1 
1 <1 k Q a l< 2 
I k 0 a l> 2 


where ko is the acoustic wavenumber as defined 
before and a is a length scale. This particular 
source definition is used here, instead of the unity 
constant amplitude, because it corresponds to a 
regular function in time domain. For unity constant 
amplitude, the time domain source would be a 
singular delta function. Since we also use a code 
based on CAA methods to solve this problem and 
the CAA code is in time domain, the regular 
source function is used to facilitate the CAA 
computation. This actually has no impact on the 
BEM calculations which are in frequency domain; 
the amplitude can be multiplied to the BEM 
solutions of unity amplitude, as a trivial step in the 
post processing of the BEM solutions. To illustrate 
the properties of the source definition (5.5), the 
source function is plotted in Figure 9, as a function 
of the non-dimensional frequency k 0 a. 

The analytical solution can be found very 
easily for an infinite plate. In this case, the 
reflection of the point source is simply its mirror 
image, which ensures the vanishing of the total 
normal velocity on the plate. Thus, the solution to 
the governing equation (5.1) and (5.4) can be 
written as 



Figure 9 Illustration of the source function 
used in the CAA computation 



Figure 10 Comparisons between analytical 
(curves), BEM (closed symbols) and CAA 
(open symbols) solutions for M = 0.2 at the 
surface location (-2, 0, 0) for a source 
located at h = 5. 


_ q(CQ) ^(jT-MCx-Xp))//? 2 | gi®) ik 0 (F- M(x-x„))//? 2 g) 

4nR* 4 ttR* 

where R is given by (2.19) and the tilde over it indicates that it is calculated with xo replaced by 
its image point 


x 0 ={0,0 ,-h}. (5.7) 

It can be seen that the solution is basically in terms of the function go defined in previous 
sections, but with the flow direction reversed. 
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The validation for the problem of sound reflection is done with comparisons between 
three methods, namely, the analytical solution (5.6), the BEM results and the CAA computations. 
In all the comparisons, the source is located above the plate at h = 5 and the mean flow Mach 
number is M = 0.2. To illustrate the comparisons, Figure 10 plots the convective derivative of the 
pressure, namely, the quantity, 


Dp 

Dt 


-ik 0 p + M ■ Wp. 


(5.8) 


This quantity is plotted here because it is the direct output of the CAA code. The results in 
Figure 10 are from the three methods at the surface location (-2, 0, 0) as a function of the non- 
dimensional frequency kOa, with the analytical solution plotted by the curves, the BEM results 
by the closed symbols and the CAA results by the open symbols. Both the real and the imaginary 
part, respectively denoted by the squares and the diamonds, are plotted in this figure and both 
show very good agreement between the three methods. 



BEM 



Re(p) | 

- 0.015 - 0.005 0.005 0.015 0.025 


Figure 12 Comparisons between analytical, 
BEM and CAA solutions for M = 0.2 on the 
plate for a source located at h = 5 for the 
real part of the acoustic pressure. 



- 0.015 - 0.005 0.005 0.015 0.025 


Figure 11 Comparisons between analytical, 
BEM and CAA solutions for M = 0.2 on a 
vertical plane through the source location 
for the real part of the acoustic pressure. 
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More comparisons are given in Figure 1 1 and Figure 1 1 for the real part of the acoustic 
pressure at the non-dimensional frequency koa = 0.93 with the mean flow Mach number of 0.2 
and the source location h = 5. Similarly to the comparisons in Figure 10, the agreement between 
the three methods is very good, except for the regions near the edges of the plate where the CAA 
solutions show discrepancies with the other two methods. This is actually expected because the 
CAA method uses absorptive boundary conditions at the computational boundaries to prevent the 
wave being reflected back to the computational domain. Thus, the waves are heavily dissipated 
in the regions near the edges, as is clearly from the lower diagrams of both Figure 11 and Figure 
1 1 . Comparisons for the imaginary part of the acoustic pressure have equally good agreement, 
though not plotted here. 

6. Sound Propagation over a Hump 

This section describes the validation of the BEM approach through a problem with 
irregular geometry, which is a hump on a plate. Without loss of generality, the hump is specified 
by a Gaussian curve of the form 




3 


2 e ~0.05(x?+x 2 2 ) 


( 6 . 1 ) 


which assumes its maximum at the origin of the coordinate system with a values of 2 and decays 
away from this maximum with a rate of 0.05. The geometry is shown in Figure 13. The 
governing equation for this problem is also given by (5.1), as in the sound reflection problem 
discussed in the previous section, and the boundary condition of zero normal velocity on the 
curved surface is now specified by 


n ■ Wp 



(6.2) 


on the surface defined by (6.1). 


Apparently, there is no analytical solution 
for such a problem so that the validation in this 
section will depend on comparisons with results 
from numerical solutions by CAA, which is also 
used in the previous section for the flat plate 
problem. The flap plate problem serves as an 
intermediate step to establish the relevance of the 
CAA code, which itself needs validation since it is 
also a numerical solution code. By first studying 
the flat plate problem that has analytical solution, 
both the CAA code and the BEM code are 
validated. A general curved surface problem is also 



Figure 13 Geometry and coordinate system 
for the problem of sound propagation over 
a hump in a uniform flow. 


20 






Re(p) | 

- 0.015 - 0.005 0.005 0.015 0.025 



Re(p) | 

- 0.015 - 0.005 0.005 0.015 0.025 


Figure 14 Comparisons between BEM and 
CAA solutions on the surface of a Guassian 
hump for the real part of acoustic pressure. 


Figure 15 Comparisons between BEM and 
CAA solutions on a vertical plane for the 
real part of the acoustic pressure. 


needed to validate the BEM approach, as in this section, because the flat plate problem, as well 
as the scattering by a sphere without mean flow, does not validate all terms in the BEM 
formulation. It is clear that the case without mean flow cannot validate any terms involving mean 
flows in the BEM formulation (2.17). It can also be seen that this same group of terms are also 
all identically zero in the flat plate case, because the mean flow is in parallel to the plate so that 
the component of the flow velocity in the direction of the plate normal, defined by M n , is 
identically zero. Thus, even though the flat plate problem has mean flow included, the terms in 
the BEM formulation are not fully validated by this canonic problem, which is why a more 
general curved surface problem is chosen to be studied in this section. 


Similarly to the other two validation cases in the two previous sections, the comparisons 
between BEM and CAA for the irregular hump problem show good agreement. To illustrate this, 
some examples are plotted in Figure 14 and Figure 15 for the real part of the acoustic pressure 
with the mean flow Mach number M = 0.2, the source location h = 5 and the non-dimensional 
frequency kOa = 0.98. The two figures are respectively for the curved surface and a vertical plane 
cutting through the source location. The good agreement is evident and again discrepancies are 
seen in the regions close to the edges of the plate, due to the use of absorptive numerical 
boundary condition for the computational domain in the CAA computations, as discussed in the 
previous section. 
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7. Far Field Formulation 

In many practical applications, the sources of the noise and the microphone locations are 
usually far apart from each other, separated by many wavelengths. In this case, the far field 
approximation can be used to simplify the computations. Thus, a form of the BEM formulation 
with far field applications is given in this section. In the notation used in previous sections, the 
far field condition means that the field point x is far away from the source point z, namely, 

|x| » |z|. (7.1) 

Thus, the usual far field limit can be applied to go(z-x) that is used in the BEM formulation and 
defined by (2.18). The far field approximation essentially neglects the variations due to z in the 
amplitude and expends the phase function in a power series. The amplitude of (2.18) is then 
approximated as 

1 1 1 

4 nR* An-]p 2 lx I 2 +(M ■ x) 2 ’ 

where the subscript °° on the quantity R denotes its value in the far field, 

Rl = V/? 2 lx I 2 +(M ■ x) 2 =1 x I -]p 2 + M 2 , 

where we have used M r to denote the component of the mean flow Mach number in the direction 
of the far field point x. It is defined by 


(7.2) 


(7.3) 


M r = M ■ x, (7.4) 

with the far field the directivity vectors defined by 


X = X / I X I . 

The phase function (2.18) can be correspondingly expanded as 

i/t 0 (iC-Mx-z ®)//? 2 

c , 


(7.5) 


(7.6) 


where the dependence of the phase function on the source coordinate z is explicitly factored out 
and the vector quantity <J> is defined by 

0 + M 

RR+m; 

Apparently, this is a function of the far field coordinate x only. 

With all these, the free space Greens function go(z-x) defined by (2.18) can now be 
written as 


So(z-x) = A 0 (x)<? 


-ik 0 z<DI jB 2 


(7.8) 
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where A 0 is an amplitude that only depends on the far field coordinate x, given by 


4)00 


l 


e >k 0 {lC-Mx )/ /} 2 


(7.9) 


The far field form of the free space Greens function allows us to factor out the large distance in 
the BEM computation to minimize numerical error. This becomes clear once the far field 
solution (7.8) is substituted into the integral solution (2.15), 


\q(z)e~^ I<!>/p2 dv- J A(z,x)g(z)ds, (7.10) 

4)00 V (z) s(j.) 

where the overhead bar on the quantity A is its normalized version given by 

A=ik 0 [2M n ~(M n + V s ■ [e'^' pl M n M s ) (7.11) 

Clearly, the right hand side of (7.10) now does not contain the large distance between the source 
and the field points. 


8. Calculation of Derivatives 

In applications of sound prediction by acoustic analogy or other acoustic theories where 
the sources are modeled as multiples, derivatives of the Greens function are needed. When BEM 
is used as a Greens function solver, the solution g is the Greens function and it is required to 
compute derivatives such as 


dx i 


and 


a 2 g 

dxfaj 


(8.1) 


Once the solutions of g are found in all spatial domains, the derivatives in (8.1) can be found in 
principle by numerical differentiation. This approach, however, has two main drawbacks. One is 
the heavy requirement for the numerical computations, for both the function g and its derivatives, 
which can be very demanding when large spatial domains are involved. It is also inefficient since 
the numerical differentiation has to be performed for all orders, even if only the highest order is 
needed in the application. The other drawback is the potential accumulation of numerical errors 
incurred in the successive numerical differentiation. These accumulated errors can easily 
overwhelm the solutions, rendering the results useless. In this section, we derive formulas for 
computing the derivatives directly from the solutions g on the body surfaces, namely, directly, 
from the outputs of the BEM matrix equation, without having to compute g itself and its 
derivatives that are not needed. 


The Boundary Element Method described in the previous sections gives the solutions of g 
on the solid surfaces with coordinates denoted by z. The solution is written as 
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(8.2) 


g(x)= f g 0 (z-x)q(z)dv- jA(z,x)g(z)ds, 

V(z) S(z) 

where the integrand of the surface integral is defined by 


A = M 0 M.g 0 + (M„M l -n I )^ + V,-(g 0 M„M,l 


(8.3) 


Here the quantities go is given by (2.18). It can be seen that the field point coordinate x only 
appears in go so that when substituting (8.2) into (8.1) to compute the derivatives, the 
differentiations with respect to x can be transferred to inside the integrations and further onto the 
function go. This leads to 


and 


jte 

dx , 


f dg o / . , f 3 A(z,x) 

J —q(z)dv- J — g(z)ds, 

V(z) 5( Z ) 



dx : dx , 

' J 


[ d2g o 
vL dz ' dz J 


q(z)dv - f 

S(z) 


d 2 A(z,x) 

dx/dxj 


g(z)ds. 


(8.4) 


(8.5) 


It can be seen that the differentiations for the free space Greens function have been transferred 
from x to z, which differ from each other by a sign. The differentiations for the integrand 
function A can also be carried out with the results 


'A = -2i k„M n (M„M, - n , ) 

oXi dZi oZjOZ/i 


a 2 g 0 _ v r^c 


dz, 


M n M s 


(8.6) 


and 


d 2 A(x,z) 

r)xr)x j 


2„ ^ 3 „ f ^2 

" Mk ~ nkJ c)zdzdz. 


r ) 2 r ) 3 

2i k 0 M n ^ - + (M n M k - n k ) , , g °. + V . • 


dZjdzj 




d 2 g o 


MM 


(8.7) 


Here all differentiations have been transferred to z, which is always the first vector in the 
combination z-x in the free space Greens function. 

With the solution g(z) given on the surfaces from the BEM solution, these results 
compute the derivatives of the Greens function. It only remains to carry out the analytical 
differentiation on the free space Greens function go, all with respect to z. With the free space 
Greens function go given by (2.18), it is easy to shown that 


dg o 

dz,- 


go F i 


where F, is the components of F, defined by (2.23), namely 


( 8 . 8 ) 
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(8.9) 


where we have introduced M, to define the components of M and the quantity A,- denotes 

, _dR* _P\z i -x i ) + M i ( z-x)-M 

' dz t R* 

From (8.8), the second order derivatives of go can be found in the form 


(8.10) 


d 2 g o 
dz/dzj 


gc 


pfj+lr 

Sz u 


By differentiating (8.9) with respect to Zj, we have 


a^ = j_ 

" ** 


f; 


V 


ifcp 1 

P 2 R 


\ 


A,.A, 

/? 


where is introduced to denote 


b u =p 2 s u +m i m j -; \ kj. 


(8.11) 


(8.12) 


(8.13) 


The above results give the first and the second order derivatives of g 0 . The process can be 
continued to find higher order derivatives. For the third order results needed in (8.7), the 
differentiation of (8.1 1) with respect to Zk leads to 




f 


d Zi dZjdz k 


= g c 


dF; 


c)F, 


FFF, +F . — j - + f. 

' 1 k ’dz k J d Zi 


+ F„ 


dF : dF : 


+ ■ 


d Zj d Zj dz k J 


(8.14) 


Here the only new term is the last one in the bracket on the right hand side and it can be found 
from (8.12) as 


dF t _[ 2 _i k 0 ' 

dZjdz k [r* p 1 

The procedure described above can be carried on to any other high order derivatives and at each 
order, there is only one new term introduced in the result. 

9. Conclusions 

In this report, we have documented the formulation, implementation and validation of a 
BEM code for computing sound propagation in uniform mean flows. The formulation follows 
the standard approach in BEM to derive an integral equation, but special attention is paid to the 
formulation of the terms that are involve the gradients of the unknown due to the mean flow. 
These terms have been reformulated in terms of the unknown itself alone so that the approach 


-r Ua 








(8.15) 
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discussed here does not have to perform numerical differentiation in the integral equation, 
avoiding potential numerical error and complication in the numerical implementation. Nor does 
it need to include the gradients as part of the unknown set in the numerical solution, and hence, 
minimizing the requirement for computation resources. The numerical implementation follows 
the approach of collocation, but sub-triangle division is used in the surface integration which 
ensures that the field point and the surface integration points never coincide. This allows the 
BEM formulation to be implemented in a very straightforward way without having to perform 
asymptotic analysis on the integrand function to avoid numerical singularity. For applications 
where the field points are far from the source points, a common situation in many practical 
applications, a far field version of the BEM formulation has been derived that explicitly factors 
out the large distance between the field points and the source points in the BEM integral 
equation, minimizing potential numerical errors. Explicit formulations have also been derived for 
computing the derivatives of the unknown function, which has applications in using BEM as a 
Greens function solver. 


Nomenclature 


A 

A 0 

C 

Fi 

Kj 

M 

N 

Q 

R* 

s 

V 

a 

C() 

8 

go 

h 

ko 

n, 

P 

q 

X 

z 


= kernel function in BEM integral 

= far field asymptotic solution for free space Greens function 

= auxiliary tensor in high order derivatives 

= factor in BEM equation to specify locations of field points 

= auxiliary tensor in high order derivatives 

= number of sides of /th element 

= mean flow Mach number 

= number of elements 

= integrated source term 

= modified radial distance 

= aggregate surfaces 

= source domain 

= radius of sphere 

= constant sound speed 

= acoustic pressure or Greens function 

= free space Greens function 

= source location above plate 

= acoustic wavenumber 

= ith component of surface normal 

= sound pressure 

= source distribution 

= field coordinate vector 

= surface coordinate vector 
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0 

po 


CO 


= coefficient of BEM matrix equation 
= auxiliary quantity involving mach number 
= emission angle in flyover plane 
= constant mean density 
= angular frequency 
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